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Abstract 

Dispersive shock waves (DSWs) in the Kadomtsev-Petviashvili (KP) equation and two dimensional 
Benjamin-Ono (2DBO) equation are considered using parabolic front initial data. Employing a 
front tracking type ansatz exactly reduces the study of DSWs in two space one time (2 + 1) 
dimensions to finding DSW solutions of (1 + 1) dimensional equations. With this ansatz, the KP and 
2DBO equations can be exactly reduced to cylindrical Korteweg-de Vries (cKdV) and cylindrical 
Benjamin-Ono (cBO) equations, respectively. Whitham modulation equations which describe DSW 
evolution in the cKdV and cBO equations are derived in general and Riemann type variables are 
introduced. DSWs obtained from the numerical solutions of the corresponding Whitham systems and 
direct numerical simulations of the cKdV and cBO equations are compared with excellent agreement 
obtained. In turn, DSWs obtained from direct numerical simulations of the KP and 2DBO equations 
are compared with the cKdV and cBO equations, again with remarkable agreement. It is concluded 
that the (2 + 1) DSW behavior along parabolic fronts can be effectively described by the DSW 
solutions of the reduced (1 + 1) dimensional equations. 

Keywords: Dispersive Shock Waves, Kadomtsev-Petviashvili Equation, Two Dimensional 
Benjamin-Ono Equation. 


1. Introduction 

In recent years the study of dispersive shock waves (DSWsjhas generated considerable interest. 
In water waves DSWs have also been termed undular bores lj Q. In fact, an early observation 
of an undular bore goes back to 1850 [sj. In plasma physics a careful observation of a DSW, 
sometimes referred to as a collisionless shock wave, was made over 40 years ago [ 3 ]- More recent 
experiments/observations of DSWs have occured in other fields, e.g. Bose-Einstein condensates 
(BEC) !H, @ and nonlinear optics Mathematically speaking the study of DSWs is difficult 

since the profile of the shock wave is highly oscillatory and the underlying shock solution does not 
converge strongly. A prototypical example of a DSW occurs in the KdV equation 

u t + uu x + e 2 u xxx = 0 (1.1) 


< 1 and inital conditions corresponding to a simple unit step (Heaviside) function. In 1974 


employing an averaging method pioneered by Whitham 0 , Gurevich and Petiavskii 


11| gave a 


detailed description of the associated DSW. About 10 years later Lax and Levermore [12| described 
the DSW rigorously via inverse scattering transform methods. Over the years there have been 
numerous important analytical studies that employ Whitham methods cf. [ll], [TTi [Til . 
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Here we study two space one time (2 + 1) dimensional equations, including the 
Kadomtsev-Petviashvili (KP) 0 and the two dimensional Benjamin-Ono (2DBO) [0] equations 
(see Eqs. m and (El), by exactly reducing these equations to the (1 + l)-dimensional 
cylindrical KdV (cKdV) and cylindrical Benjamin-Ono (cBO) equations (see Eqs. (12.181) and (12.191) 1 
respectively. 

We analyze the cKdV/cBO equations via Whitham theory and derive the general Whitham 
modulation equations; these equations are transformed into simpler form by introducing appropriate 
Riemann type variables. These Whitham equations in Riemann variables are not in diagonal form. 
We remark that in the cKdV case diagonal form may be obtained using the integrability of cKdV 
on the other hand, neither 2DBO nor its reduction, the cBO equation is known to be 
integrable. 

We study the DSWs in the cKdV and cBO equations numerically and describe their differences 
from the DSWs in the classical KdV and BO equations. Indeed the DSWs in the former are found 
to decay slowly in time whereas those in the latter do not exhibit such temporal decay. We find 
that direct numerical simulations of the Whitham modulation equations agree well with those of the 
cKdV and cBO equations. 

We then compare these (1 + 1) dimensional DSW structures to direct numerical simulations of 
the (2 + 1) dimensional KP and 2DBO equations. Our comparisons between 1 + 1 numerics/theory 
and 2 + 1 numerics exhibit excellent results. In general the DSW weakens across the parabolic front 
as time increases. We also note that the numerical simulations of 1 + 1 Whitham theory which 
removes the fast variation, are much faster than the 1 + 1 cKdV/cBO equations which in turn, are 
orders of magnitude faster than the 2 + 1 equations. 

Over the years there have been many numerical studies and calculations associated with the KP 
equation cf. [2214261 ]. 

Our interest is to study DSW systems which have a discontinuity across a parabolic front; this is 
analogous to the well known Riemann or shock tube problem in classical shock waves. We find that 
indeed there are DSWs generated across the shock front. To our knowledge this is the first time the 
nondecaying 2 + 1 analogue of a Riemann problem for KP and 2DBO is analyzed in detail. 

The reduction discussed here, which we term parabolic front tracking, was employed as 
self-similar reductions [ 27 ], 0] in the analysis of the Khokhlov-Zabolatskaya (KZ) equation [0] 
(see also 0). Indeed the KP/2DBO equation reduces to the KZ equation in the limit of zero 
dispersion. When viscosity is added to the KZ equation the relevant shock waves are strongly 
convergent. Our study of the KP/2DBO DSWs requires critical use of Whitham modulation theory, 
which is necessary due to the weak convergence of the DSWs. 

This paper is organized as follows. In Section 2 we reduce the KP and 2DBO equations to the 
cKdV and cBO equations along a parabolic front. In Section 3 we employ perturbation theory 311 
to find the conservation laws associated with Whitham theory for the KdV and cKdV equations. We 
then transform the Whitham modulation equations employing Riemann-type variables; the resulting 
Whitham system is not immediately diagonalizable. We solve the 1+1 Whitham system associated 
with the KdV and cKdV equations numerically and reconstruct the DSW solutions of KdV and 
cKdV. We then compare these results with direct numerical simulations of KdV and cKdV and 
show that, apart from an unimportant phase they are in excellent agreement. We also note that the 
Whitham equations for cKdV exhibit a small discontinuity. This discontinuity would be resolved by 
taking into account higher order terms (see 0), but doing so is outside the scope of this paper. In 
Section 4 the BO and cBO equations are analyzed in the same way as KdV and cKdV are analyzed 
in Section 3. In Section 5 we compare the 1+1 results for cKdV/cBO and the 2+1 results for 
KP/2DBO by direct numerical simulations. After accounting for an unimportant mean term we 
again find excellent agreement; animations are also included as part of our 2 + 1 description. We 
conclude in Section 6. 
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2. Reduction of KP, 2DBO equations to cKdV, cBO equations 


In this section, we examine DSW propagation associated with two different (2 + 1) dimensional 
nonlinear partial differential equations (PDEs). One is the Kadomtsev-Petviashvili (KP) equation 


(ut T UU X + € Uxxx) x T Au^y — 0 


( 2 . 1 ) 


where e, A are constant. This equation was first derived by Kadomtsev-Petviashvili 17] in the context 
of plasma physics; subsequently it was derived in water waves [33j where it describes the evolution 
of weakly nonlinear two dimensional long water waves of small amplitude. When |e| <C 1 we have 
weak dispersion. According to the sign of A, Eq. (12.11) is usually termed KP-I (—) or KP-II (+), 
respectively. KP-I describes the dynamics when the surface tension of the water is strong and KP-II 
describes the dynamics with weak surface tension. The other equation we study is 


(ut + UUx + PP {^xx')') x T A Uyy — 0 

where Pu(x) denotes the Hilbert transform: 


( 2 . 2 ) 


Hu(x) = -V f dx' (2.3) 

7 T Jx' - x 

and V denotes the Cauchy principal value. We refer to Eq. (12.21) as the 2DBO (Two Dimensional 
Benjamin-Ono) equation; it is a two-dimensional extension of the classical BO equation and describes 
weakly nonlinear long internal waves in fluids of great depth flil ]. 

The goal in this paper is to enhance understanding of DSWs in multidimensional systems. A 
general form for these two equations is 


(Ut + UUx + -^(^))x A Uyy 0, (2.4) 

when Fi(u) = e 2 u xxx Eq. (12.41) is the KP equation and when ^(it) = PH(u xx ) it is the 2DBO 
equation. Also, we are interested in a class of initial conditions for Eq. (IQ) which are almost 
step-like initial data, for example 

1 + y tanh ( K 


x +2 P (2b°) 


(2.5) 


u(x,y, 0) = - 


where y — ±1 ,K are real constants. 

The above front (12.51) is a regularized two dimensional extension of the Riemann-type initial 
condition 


u(v, 0) 


1 , r] < 0 ; 

0 , r] > 0 , 


( 2 . 6 ) 


where r/ = x + \P(y , 0). For the special choice of a parabolic front 


P(y, 0) = q/ 2 , 


(2.7) 


where c is a real constant, graphs of the initial conditions with y, = 1 (y, = —1) are given in Fig. [I] 
Note that these initial data increase (decrease) in x across the parabolic front. 

We assume that solutions of Eq. ( 1 P 1 ) satisfy the following ansatz: 


u = f{x + P(y,t)/2,t,y); 


( 2 . 8 ) 


the front is then described by x + P(y , t )/2 = constant. We substitute the ansatz (12.81) into Eq. (12.41) 
and find 


Ptfri + ft + / ft] + Pi (fiv)) ) 


4 {Py) fm + ^Pyyfv + Pyfvv + fyy ) — 0 (2-9) 
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Figure 1: Contour plots of initial conditions for P(y, 0) = cy 2 where c = 0.1, K = 10, and (a) y = 1 
(non-decreasing in ie); (b) y = — 1 (non-increasing in x). 


where r] = x + P(y, t)/ 2. In the above equation u satisfies the following boundary conditions at the 
infinities: for non-decreasing type initial conditions 

u —> 0 as r\ — > — oo and u — >• R(t) as p — > oo (2-10) 

and for non-increasing type initial conditions 

u — > R(t ) as 77 — ► —00 and u — > 0 as rj —> 00. (2.11) 

The function R(t) is chosen appropriately; the forms of R{t) are given later in this section. 

Using these boundary conditions and assuming that P yy is independent of y, the coefficient of 
the term f v vanishes in Eq. (12.911 and thus / can be taken to be independent of y. Then we obtain 
the following system of equations 

A 9 

Pt + — (P y ) =0, (front shape equation : FS ) (2.12a) 


ft + f fj ] + 77 Pyyf + Fi (f(v)) — 0- 


(2.12b) 


The FS equation can be transformed into the inviscid Burgers equation, also sometimes called the 
Hopf Equation 341: 

Vt + A vv y = 0, (2-13) 


by taking v = P y . This equation is exactly solvable by the method of characteristics. For the 
parabolic front initial condition (12.71) . we get the following initial condition for Eq. (12.131) : 


v(y, 0) = 2 cy. 


The solution of the initial value problem (IVP) (12.131) and (12.141) is found as 


v(y,t) 


2cy 

1 T 2 cAt 


Thus the front shape P(y,t) is given by 


P(y,t) 


cy 


2 


1 + 2cXt ' 


(2.14) 


(2-15) 


(2.16) 
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Consistent with our original assumption, the front curvature P yy is indeed independent of y for all 
time. 

At fixed time t, the level curves of the solution rj = x + p ^ i ' > = ?yg, rjo constant, are along the 
following curves in the (x, y)-plane: 


x = — 


~ 9 

cy 


2(1 + 2 Act) 


+ Vo = C{t)y 2 + ?7o, C[t) = — 


2(1 + 2Acf)' 


(2.17) 


This shows that when cA > 0 the curvature of the initial parabolic front decreases in the positive 
t direction, or equivalently the parabolic front ‘flattens’ when t increases as indicated by C(t). 
However, the curvature blows up in the negative t direction at a critical time t c = —1/(2Ac). In the 
following we will take cA > 0 and only be concerned with t > 0 to avoid blowup. Evidently this 
reduction is not valid for all time, but this is common with self-similar reductions as they frequently 
describe asymptotic phenomena. 

In summary, we have shown that for the parabolic front initial condition (12.71) . by using the 
ansatz (12.81) . the (2 + 1) dimensional PDE (12.41) can be exactly reduced to a (1 +1) dimensional PDE 
(12.12bl) with variable coefficients.This 1 + 1 dimensional PDE is either the cylindrical KdV (cKdV) 
equation 

ft + f f v + ^ + 2\ct^ + = 0 ( 2 - 18 ) 

or the cylindrical BO (cBO) equation 

ft + ffv+ 1 + A 2 C Aa f + en (M = 0 (2-!9) 

depending on the choice of F,j (u) in Eq. (E3). We reiterate that the solution (12.161) shows that the 
term P yy in (12.12b () is independent of y. In the next two sections, we examine the DSW solutions 
of Eqs. (12.181) and (12.191) with (non-increasing) initial data such as Eq. (12.61) . 

Later we will denote to = so the term 1 ^ X 2 c Xdt becomes r 2 t+t 0 ) ~ Also, we will consider only 
A = 1 (i.e,. KPII, 2DBOII); the other sign (KPI, 2DBOI) can be obtained by changing c to — c, i.e. 
changing the direction of the parabolic front. 

In order to determine the boundary conditions associated with Eq. (12.18112.191) at infinity, we 
neglect rj dependent terms and then solve the remaining ODE with the corresponding initial condition 
(ESI). The solution of this ODE with the initial condition R( 0) = 1 determines the function R(t) in 
the boundary conditions (12.101) and (12.111) as 


m 


i 

\/l + 2 ct 


( 2 . 20 ) 


We note that when the parabolic front initial condition m is not satisfied, the term P yy 
in (I2.12bl) is no longer independent of y, and so Eq. (12.12bl) cannot describe the front evolution 
self-consistently. In this case, taking the FS equation (12. 12al) to be the same, Eq. (12. 12bl) then 
changes to 

ft + TjPyyf + / fv + + HPyfvv + fyy) = ( 2 - 21 ) 

This new equation is a (2 + 1) PDE which at this point is more complicated than the original 
equation (1241) . Hence the ansatz (12.81) for general initial conditions does not lead to cKdV or cBO 
since there are additional terms which depend on f yi f yy - However for a class of initial conditions 
whose evolution in time satisfies the condition 


|L(/;P)|«1 (2.22) 

where L(/; P) = Pyfmj+fyy: the system (12.121) still approximately describes solutions of the equation 
(12.41) perhaps for a finite time and/or a subdomain in space. The condition (12.221) means the 
additional terms in Eq. (12.211) can be neglected and Eq. (12.211 ) turns into Eq. (I2.12bl) (with A = 1). 
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3. Dispersive shock waves in the KdV and cKdV equations 


In this section we will investigate the DSW solutions associated with both the KdV and cKdV 
equations via Whitham modulation theory. We will compare Whitham theory and direct numerical 
simulations and show that they agree well. Since the KdV equation may be viewed as a special case 
of the cKdV equation (12.1811 when c = 0, we will develop the Whitham approach only for cKdV. 

From Whitham modulation theory we find three conservation laws. Then, we transform these 
three conservation laws into a system of quasilinear first order PDEs by using convenient Riemann 
type variables; this was first introduced/derived by Whitham for the KdV equation |Toj]. For KdV 
this system can be diagonalized and solved exactly. The system for the cKdV equation cannot be 
immediately diagonalized. Numerically we show that the system has solutions which demonstrate 
the DSW structure of cKdV, unlike KdV, decays in time. 

We will use a method of multiple scales originally employed by Luke [31] in the study of Whitham 
type systems associated with a nonlinear Klein-Gordon equation. For the cKdV equation the leading 
order equation has a Jacobian elliptic function solution, i.e. the cnoidal wave solution, where the 
parameters are slowly varying; there are three independent parameters. The leading order problem 
introduces the rapidly varying phase which requires a compatibility condition which is often termed 
conservation of waves. The next order in the perturbation method has two secularity conditions; 
these together with conservation of waves give three conservation laws. 


3.1. Whitham modulation equations for KdV/cKdV - conservation laws 

In what follows we develop the slowly varying Whitham modulation equations. We assume 
/ = /( 6 *, 77 , e) where 9 is rapidly varying and defined from 


k(r],t) ^ = u(v,t) = kV 


(3.1) 


where k, to and V are the wave number, frequency and phase velocity, respectively. This definition 
gives us the compatibility condition (9 v ) t = (0*) (conservation of waves) as 


k t + (kV) v = 0. 

This is the first conservation law. 

With these rapid and slow variables we transform Eq. (12.181) using 

8 k d d d oj d d 

dr, ^ e 89 ^ dr,’ dt * e 89 dt 


(3.2) 


(3.3) 


we have 


or 


,, oj d d 2 /k d d 3 .fed d 

[(_ 7 dO + di )+C i ~ed9 + fr ] ) ]f + f{ ld0 + fr) )f 


Ac 


1 + 2A ct 


f = 0 


1 ( df df 
— oj-^ + kf- 


d9 


89 


V vv 


df 

89 


k 503 ) + 

d 2 f , , d 3 f 


' 89dr) 89dr) 2 


df_ j.df_ 
dt 1 dr) 

I 2/c ^ f ) 

d0 2 8r)) 


3 kk. 


3 k 


d 2 f 

' n 89 2 

e 2 ^ = o. 

8r, 3 


2 d 3 f 


Ac 


89 2 dr, 1 + 2A ct 


(3.4) 


Then we expand / in powers of e as 


/ {&, 0, t) = f 0 (0, r), t) + e/i (0, 77 , t) + .... (3.5) 

Grouping the terms in like powers of e gives leading and higher order perturbation equations; we 
only consider the first two orders here. The O (■)■) equation is 

-w/ 0,0 + kf 0 fo,e + k 3 f 0 ,888 = 0; (3.6) 
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the 0(1) equation is 


•C/l = —ufl,9 + k (fofl)g + k 3 fi,ggg — G, 


(3.7) 


where 


G = — I /o,t + /o/o,i7 + Skk^fo^gs + 3fc 2 /o,eer) + 


fo 


2 1 + to 


We can proceed to higher order terms, but doing so is outside the scope of this paper. 
The solution of (13.61) is 

fo (0 , rj,t)=a ( 77 , t) + b (r,, t) cn 2 [2 (9 - 9 0 ) K, m ( 77 , t)] 
where K = K (m ( 77 , t)) is the complete elliptic integral of the first kind and 

b b 


k z = 


48m 2 K 2 ’ 


a = V 


_ 

3m 2 3 ’ 


(3.8) 


(3.9) 


(3.10) 


and recall that V = f. For the purposes of this paper we consider do to be a constant. At this point 
we have three independent parameters: b , V. m. 

We rewrite the conservation law (13.21) by using the above formulae as 


9 1 

< 1 ,m 

+ 9 1 

( v J 

dt 1 

l 4 V3K \m 2 

+ dr, 1 

{ aVsk V 


= 0 . 


(3.11) 


When we use the solution m in (EZD, secular terms can occur, i.e. terms that grow arbitrarily 
large with respect to 9. Let w denote solutions of the adjoint problem to Cu = 0, i.e., 

(3.12) 


C A w = 0, C A = ujdg — kf 0 dg — k 3 c 
To eliminate the secular terms, we use the following relation obtained from Eq. (1X71) 


[ [wCfi — f\C A w\dO = f wGdO. 
Jo Jo 


(3.13) 


The adjoint problem (13.121) has two linearly independent solutions w = 1 and w = fo, the latter 
following from Eq. (13.61) . We substitute these into Eq. (13.131) . enforce the periodicity of fo ( 6 , 77 , t) in 
9 and obtain the following secularity conditions 


Gd9 = 0, and / foGdd = 0. 


(3.14) 


Using (cf. 


for i = 1 , 2 ,... and j = 1 , 3, ..., and 


d9 i 


d9 = 0, 


''O'- 


/o 


f fofo,eed9 = - [ flgdd, 
Jo Jo 


we get from the first secularity condition in (13.141) 

fodO 




dt 


1 


2 1 + to Jo 


fod9 = 0 


and the second secularity condition in (13.141) 


d_ 

dt 


(3.15) 


(3.16) 


(3.17) 


rlfSde+ i{U - “ 2 1', fl ' de ) + *rrJ! f ° de = a 
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(3.18) 



















We can use Eq. ED to rewrite the conservation laws (|3.17l) and (13.181) in terms of to, V and 6 /to 2 . 
From the properties of the elliptic functions 


we find 


rl b 

fodO = V + 


/ o 


3to 2 V K 


3 E 


to — 2 


l fSM = 1/2 + 2V 3^ (¥+ ml - 2 ) + (i) K - +'). 

r +m2 ~ 2 ) + f (J?) (™ 4 -"> 2+1 ) 


f*dd = V+V 


(3.19) 


1 

5 V to 2 


3 r 


E 1 

— (to 4 — to 2 + l) H-(5 to 6 — 21to 4 + 33m 2 — 22) 

K \ ) 27 v ' 




'0 


— (to 4 — to 2 + l) — (to 4 — 3to 2 + 2) 


where K = K[m) and E = E(m) are the complete elliptic integrals of the first and second kind. 
Using the formulae (13.191) in (13.171) and (13.181) the following conservation laws are obtained 


d_ 

at 


V 


+ to 2 — 2 


ld_ 

2 drj 


3 E 


3m 2 \ K 


o b f 3E 9 

V 2 + 2V—r — + to 2 - 2 


3to 2 \ K 


) + ( i ) (™ 4 -’" 2 + 1 ) 


(3.20) 


+ 


1 


and 


d_ 

at 


2 t + to 

b 


U + 


3 E 


3to 2 \ K 


+ to 2 — 2 


= 0 


U 2 + 2U 


3 E 


3to 2 \ K 


+ TO 2 - 2 ) + (to 4 - TO 2 + 1) 


d_ 

+ ~bh) 


2U 3 2U 2 6 /3 E 


3 


2F 


ra z 


4 _2 

m 


1 ) 


+s(^) ( 2 ™ 2 -1) k+1) k - 2) 

v ' 2 + 2 i/ i(¥ +m 2 ’ 2 ) + (i) (™ 4 -“ 2 + i ) 


(3.21) 


= 0. 


2t + to 

Equations (13.111) . (13.201) and (13.211) are the three conservation laws which determine 6, to and V. 

3.2. Whitham modulation equations for KdV/cKdV - Riemann type variables 

We transform these conservation laws by making a suitable change of variables. In particular we 
parametrize 6, to, V in terms of the following Riemann type variables r\, r 2 , 


6 = 2(r 2 -ri), to 2 = —-—, V = - (ri + r 2 + r 3 ), n < r 2 < r 3 

r 3 -r 1 3 

and it follows from Eq. (13.101) that a = r 3 — r 2 + r±. 

This leads to an equation of the form 


(3.22) 


6 f) r . . 

+ Ci t j) = 0 , * — 1 , 2 , 3 , 


j=i 


(3.23) 




















































where Ajj , Bij are functions of ri, i = 1,2,3 and Cij is a function of n, i = 1,2,3 and t. Multiplying 
by the inverse of the A matrix we can simplify Eq. (13.2311 into the following quasilinear PDE system 


where 


Vi 


dn dri hi (n, r 2 , r 3 ) 

- M +v i (r l ,n,r 3 )-^+ 2( + (o =0, 


1 . 2 K(m) 

o ( r i + r 2 + r 3 ) - - r 2 - rq) —— -——, 

3 3 A (in) — E(m) 


*=1,2,3, 


(3.24) 


1 2 

v 2 = - (ri + r 2 + r 3 ) - - (r 2 - rq) (l - m 2 ) 

1 . , 2 . , K(to) 

W + ^2 + r 3 ) + - (r 3 - rq) 

3 3 A(m) 


K (to) 


E(m) — (1 — to 2 ) K(m) ’ 


/*! = 
/12 = 
h 3 = 


(5E(m) — 3 K(m))ri — ( E(m ) + K(m))r 2 + (K(m) — E(m))r 3 
2>(E(m) — K(m )) 

-E(to) (r 3 - ri) (n - 5r 2 + r 3 ) - A'( TO ) ( r 2 - r 3 ) (ri + 3r 2 - r 3 ) 
3 [E(m)r i — K(m)r 2 + ( K(m ) — £?(m))r 3 ] 

(2 K(m) — E(m))r 2 + (5 E(m) — 2K(m))r 3 — E(m)ri 
3E(m) 


(3.25) 


The r,’s are called Riemann variables. From Eq. (13.22|) the solution of the leading order problem 
is reconstructed from the the r^’s as 


/o (0, f) = ri - r 2 + r 3 + 2 (r 2 - n) cn 2 [2iC (6» - 6 0 ), to] . 
The rapid phase 9 is determined by integrating S3 


Hn,t)=r k -AA dz t' 

J—00 ^ j 0 


Hv,t ) 


dt 


(3.26) 


(3.27) 


We also note that there is a free constant 9 0 in Eq. (13.261) which we determine by comparison with 
direct numerical simulations. 

The initial values of the Riemann variables of the reduced diagonal Whitham system (13.241) are 
given below (see Fig. 


n( 77 , 0 ) = 0 , 7-2(77,0) = 


77 < 0; 
77 > 0, 


7 - 3 ( 77 , 0 ) = 1 . 


(3.28) 


In the absence of cylindrical terms, i.e. hi = 0, Eq. (13.241) reduces to a diagonal system that agrees 
with the well known diagonal Whitham system for the KdV equation first derived by Whitham |lp] 
(see also llj). For the KdV equation, the solution of the reduced Whitham system provides the 
dispersive regularization for the initial data (12.61) . This regularization can written in terms of a 
similarity variable £ = 77 /t and the system (13.241) reduces to 


(*>2(7-2) -£>2(0 = 0 


(3.29) 


from which we get the self similar solution v 2 (r 2 ) — £, = 0 or by inversion r 2 = r 2 (^). Thus the system 
(13.241) with hi = 0 admits an exact rarefaction wave solution in terms of the self-similar variable 
£ = 77 /^. This rarefaction wave solution leads to the DSW solution for the classical KdV equation 
(for additonal details cf. i)- 

However, the system (13.241) for cKdV is not diagonal and this property makes finding analytical 
solutions more difficult. The solutions of nondiagonal quasi-linear systems obtained via Whitham 
modulation theory were investi gat ed in certain cases. In [l3j], stationary solutions were obtained in 
the KdV-Burgers equation. In 14|, the leading and trailing edge structures of DSW solutions were 
investigated in a variable coefficient KdV (vKdV) equation. 
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Figure 2 : Initial values (13.281) for Riemann variables n, ri and r 3 . 



3.3. Comparison between numerical simulations of Whitham modulation equations and direct 
numerical simulations of KdV/cKdV 

Our approach is to study the cKdV equation, and the KdV equation as a special case, by using 
numerical methods to solve the (nondiagonal, in general) Whitham modulation equations (13.241) . 
We compare the results to direct numerical simulations of the original 1 + 1 cKdV equation. This 
allows us to understand the underlying structure of the DSWs in the cKdV equation. Indeed 
we find very good agreement between the numerical solutions of the Whitham equations and direct 
numerical simulations of the cKdV equation. We conclude that Whitham modulation theory provides 
a good approximation of the DSWs in the cKdV equation. The advantage of computing with the 
Whitham system is that it gives the structure of the DSWs in terms of 0(1) coefficients, whereas 
for direct numerical simulations one has small coefficients (due to e 2 1 ) which in turn requires 
more sophistication to solve and longer computing times. 

First we find numerical solutions of the diagonal reduced Whitham system associated with the 
KdV equation and the nondiagonal Whitham system (13.241) for the cKdV equation; the initial values 
of Riemann variables are given by Eq. (13.281) . 

The boundary conditions for Eq. (13.241) must be determined before numerical computations can 
proceed. For the KdV equation and its associated Whitham system, the boundary conditions remain 
constant at both ends of the domain. However, the boundary conditions change in time for both the 
cKdV equation and its associated Whitham system (13.241) . The boundary conditions for the cKdV 
equation are the same as in Eq. (12.111) . For the Whitham system (13.241) . we get the corresponding 
boundary conditions by numerically solving the reduced ODE system obtained from Eq. (13.241) by 
neglecting the spatial variable 77 . This ODE system is solved with the initial conditions (13.281) at both 
ends separately by using the ode45 solver of MATLAB®. The evolution of the Riemann variables 
for the cKdV equation at the boundaries is given in Fig. [3] We see that these Riemann variables at 
the boundaries decay in time. 

For the computation of the Whitham system (13.241) including boundary conditions, we use a first 
order hyperbolic PDE solver based on MATLAB® by Shampine 37j and choose a two-step variant 


of the Lax-Wendroff method with a nonlinear filter 


In the numerical solutions of Whitham 


systems, we use N = 2 12 points for the spatial domain [—30,30] with the time step being 0.9 times 
the spatial step. The resulting numerical solutions are given in Fig. [4] for both KdV and cKdV. 

Interestingly, in the cKdV case the r$ component of the solution of the Whitham system exhibits 
a small shock-like front in front of the DSW. Direct numerical simulations indicate that this behavior 
is not significant. In fact adding higher order terms to the Whitham system is expected to regularize 
the solutions (see [321]'). 

We can reconstruct the corresponding DSWs at any time (e.g. t = 7.5) for both KdV and cKdV 
from the Riemann variables rfs using Eqs. (13.261) and (13.271) . These are also plotted and compared 
with direct numerical simulations of KdV/cKdV (discussed below) in Fig. [5^ and Fig. [5)). In Fig. [SJo 
we have chosen the arbitrary phase @0 i n Eq. (13.261) appropriately to agree with the direct numerical 
simulations. 
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Figure 3: Evolution of Riemann variables for cKdV case which are obtained by numerical solution of reduced ODE 
system (13.241) (a) at the left boundary, (b) at the right boundary. 




Figure 4: Riemann variables at t=7.5 which are found by numerical solutions of Whitham system (13.241) (a) 
for the KdV eq., (b) for the cKdV eq. Here, we take to = 10. 


For the direct numerical simulations of KdV/cKdV, we use a numerical procedure which is useful 
for problems with fixed boundary conditions. Since the left boundary condition R(t) for cKdV is a 
function of t we first transform (12.181) by 


/ = R{t)</> 

to the following variable coefficient KdV (vKdV) equation 

R(t^) 4*4*7] “f C 4rir}r) — 0; 


(3.30) 

(3.31) 


where we recall from Eq. (12.201) that R(t) = J^_ to with to = 
boundary condition fixed at <(>_ = 1 , while the right boundary 
as in the original cKdV equation. In order to solve Eq. (13.311) 
differentiate with respect to 77 and define 4*v = z to get 


= i. Equation (13.311) has the left 
condition (f> + = 0 stays the same 


numerically (see also 3§-41|) 


we 


z t + R(t) ( z<f)) rj + e 2 z vvv = 0. 
Transforming to Fourier space gives 


(3.32) 


St = Lz + R(t )N (z, t ) 

where ~z = T(z) is the Fourier transform of z, Lz = ie 2 k 3 z and 


(3.33) 


N (£, t) = — ikJ- 


(z) dr] 


T 


-1 



(3.34) 
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Figure 5: Numerical and asymptotic solutions of KdV and cKdV eqs. at t=7.5 with the initial data COO): 
(a)for KdV eq., (b) for c.KdV eq. Here, we take to = 10 and e 2 = 1CP 3 . 


The only difference from the classical KdV case is that for vKdV the nonlinear term N has a time 
dependent coefficient (see (39[). To solve the above ODE system in Fourier space we use a modified 
version of the exponential-time-differencing fourth-order Runge-Kutta (ETDRK4) method mm. 
For the required spectral accuracy of the ETDRK4 method, the initial condition for z must be 
smooth and periodic. However, the step initial condition (12.61) for u or equivalently / leads to 
^( 77 , 0) = — 15 ( 77 ), where <5 represents the Dirac delta function. Therefore we regularize this initial 
condition with the analytic function fiiol 


z(t 7 , 0 ) = ^-sech 2 [kr^j , 


(3.35) 


where K > 0 is large. Thus Eq. (13.321) can be solved numerically via Eqs. (13.33113. 3T1) on a finite 
spatial domain [—L,L\, where T represents the discrete Fourier transform and the integration limit 
—00 in (13.341) is replaced by — L. 

For the ETDRK4 method we take the number of Fourier modes in space to be N = 2 12 , 
the domain size to be L = 30, and the time step to be 10 -4 . The parameters are chosen to be 
cT 1 = to = 10, e 2 = 10 ~ 3 and K = 10. The numerical results for the KdV equation (with c = 0) and 
the cKdV equation are given in Fig. [5jr and Fig. [5}o at t = 7.5. To provide another view, we also 
include space-time plots of the direct numerical solutions of KdV and cKdV eqs. are given in Fig| 6 ] 
For both KdV and cKdV, the structure of the DSWs, its leading edge amplitudes and the 
wavelength of the oscillations predicted by the asymptotic solutions agree well with the numerical 
solutions. However in both cases, the position of the leading edges are slightly different, i.e. there 
is a small phase shift; this phase shift is larger towards the rear end of the DSW. Indeed one needs 
to proceed to higher order terms in the asymptotic expansion to achieve better results for the phase 
[32j and to smooth any shock-like discontinuities. This topic is outside the scope of this paper but 
we will return to it in a future communication. In our comparisons we have chosen the phase shift 
in the reconstructed leading order solution (13.261) to adjust the first maxima to agree with direct 
numerical simulations. This phase shift is not fixed due to the asymptotic nature of the initial value 
problem. 

For the KdV equation, the amplitude of the trailing edge remains fixed at 1, while the amplitude 
of the leading edge approaches 2. The leading edge position is 77 + = 4.946 at t = 7.5, so the average 
front speed of the DSW is calculated to be V+ = 0.659. This speed agrees almost perfectly with the 
phase speed V+ = 2/3 of the soliton solution of the KdV equation with amplitude 2: 


fiVi t) = 2sech 2 


. -\/6e 2 


V-Vo--t 


(3.36) 


For the cKdV equation, the amplitude decays in time for the trailing edge. Accordingly the 
amplitude of the leading edge also decays in time, and so does the speed of the leading edge. Indeed 
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Figure 6: Space-time plot of the direct numerical solutions between t = 0 and t = 20 (a)for KdV eq., (b) for 
cKdV eq. Here, we take to = 10 and e 2 = 10 -3 . 


in this case the leading edge position is rj + = 3.549 at t = 7.5, so the average front speed of the 
DSW is approximately V+ = 0.47. This is smaller than the average front speed in the KdV case. 

Finally we remark that the numerical solution of the Whitham KdV equations in Fig. 2^, suggests 
that there is a self-similar solution of the reduced Whitham system (13.24|l when hi = 0; indeed this 
is true and well-known. We note also that in general the Whitham system (13.241) associated with the 

cKdV equation has a similarity solution of the form ri(rj,t) = r,; (, i = 1, 2,3 where £ = 77 /(2 t+i 0 ); 


the similarity equations are given by 


drj 

dl 


(vi (n , r 2 , r 3 ) - 2gJ + hi (n, r 2 , r 3 ) 


= 0, 


i = 1,2,3. 


(3.37) 


However this similarity system is unlikely to be uniformly valid on the whole domain. 


4. Dispersive Shock Waves in the BO and cBO Equations 

In this section we will study DSWs associated with the BO and cBO equations. We will follow 
the method described in the above section for the KdV and cKdV equations. 


f.l. Whitham modulation equations for BO/cBO - conservation laws 

We begin with Eq. (I2.19[) and introduce the transformation of variables (13.31) into fast and slow 
coordinates. This yields the equation 


r.Lod d 






Ac 


1 + 2 Act 


/ = 0 


or 


~ u aS + kf m 

df , M 


k 2 n 


dt 


+f df ] +n \ k 


dj_ 

'86 


de 2 
+ 2 k 


d 2 f 

09 drj 


+ 


Ac 


1 -|- 2A ct 


eU 


Orf 


(4.1) 


= 0. 


We then expand f = fo + e/i + ... with the equations of the first two orders given at 0 ( 7 ) by 

—ufo,e + kfofo,e + k 2 H {fo : ee) = 0; (4.2) 
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and at 0{ 1) by the linear equation 


£ff/i = —w/i t g + k (u 0 /i)g + k :2 'H {f\.ee) = G 


where 


G=- 


fo,t + fofo,v + 'H (k v fo,8 + 2kfo t g v ) + 


fo 

2 t + to 


The solution of Eq. (Oil is 


fo (<?, V, t) 


4fc 2 


vM 2 + 4 k 2 — A cos {9 — 6q) 


(4.3) 

(4.4) 


(4.5) 


where do is constant. 

Actually, Eq. m is the periodic wave solution of the classical BO equation which was first 
obtained by Benjamin [42j]. Here A = | {fo,max — fo,min ) is the amplitude of the wave (cf. [lij) and 
the phase velocity of the wave is given by 

V= i\/A 2 +4fc 2 + /3. (4.6) 

In Eq. (14.51) . k, A, [5 and V are functions of slow variables ry and t. In what follows we get the 
general modulation equations for the cBO equation in terms of the three variables k, V and /?; note 
from Eq. Oil A can be written in terms of these variables. 

As with KdV and cKdV the conservation of waves (13.21) is a necessary compatibility condition 
for the cBO equation as well. This is the first conservation law. The other two conservation laws are 
obtained by eliminating the secular terms at the right-hand side of Eq. (14.31) . In a similar manner 
to KdV/cKdV above, let w denote solutions of the adjoint problem to Cru = 0, i.e., 

Cjjw = 0, Cjj= udg - kfod e - k 2 U {dee) (4.7) 


where we used the anti-symmetry of the Hilbert transform: (Hu,v) = (u, — T-Lv),f) being the 
standard inner product. In order to eliminate secular terms, we use the following relation that 
follows from (14.31) 

/* 27 t r>2ir 


/ [wCnfi - fiCiw]dO 


wGdd. 


(4.8) 




Jo 


We put w = 1 and w = fo (obtained from Eq. (14.21) 1 into Eq. (14.81) . enforce the periodicity of 
fo ( 9 , T], t) in 9 and obtain the secularity conditions respectively as 


Gd9 = 0, and / foGdd = 0. 


(4.9) 


Using following identity, 


7 n (w ]d ° =o ’ 


we get from the first secularity condition in Eq. (SU) 


d_ 

di. 


fod9 


d (\ 


dr) \2 


fldO ) + 


1 


2 1 + to Jo 


f 0 d9 = 0 


(4.10) 


(4.11) 


and the second secularity condition in Eq. (14.91) 


d f 27r d /2 r 271 \ C 2n 2 

ml f « m+ m,{zl f ° m ) +2 l )<ie+ — l = o. (4.i2) 
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From the properties of the Hilbert transform 43]], we have 

/ f Q dd = 2n (/3 + 2k ), 

Jo 

p2n 

/ f^dO = 2tt (f3 2 + 4Vk) , 

Jo 

/ f$dd = 2t r [/3 3 + 6/3 2 fc + 12fc/3 {V-/3) + k (3 H 2 + 8fc 2 )] , 

Jo 

p2ir 

/ foU (k v fo : e + 2kfoj v ) dO = - [kA 2 ] . 

Jo 

Substituting the definition (14.61) and identities (14.131) into Eqs. (14.111) and (14.121) we can simplify 
them to find the following conservation laws 


(4.13) 


Pt + PPr) 


2k+ /3 

2t + to 


= 0 


and 


V t + W v + kk v + 


2V — [3 


= 0. 


(4.14) 


(4.15) 


v 71 ' 2 t + t 0 

Equations (13.21) . (14.141) and (14.151) are the three conservation laws for the three variables k. V and 

P- 

4-2. Whitham modulation equations for BO/cBO - Riemann type variables 
It is convenient to introduce Riemann type variables a, b, c [16 ] 


k = b — a, V = b + a, f3 = 2 c 
and write the leading order solution fo in terms of a, b , c 


fo (0, r ), t) = 


2 (b — a) 


(b + a — 2c) — 2 y/(a — c)(b — c) cos {9 — 6q) 


+ 2c 


(4.16) 


(4.17) 


where the phase 9 is determined by Eq. (13.271) and 9q is at this stage an arbitrary constant. 

This transformation to Riemann type variables simplifies the conservation laws. We can write 
the quasilinear PDE system for Riemann variables a, b and c in the following form (recall fy = 1/c) 


„ a + b — c 

at + laciri H———j—-— — 0, 


2 bb n 


Ct + 2cc„ + 


2 1 + to 
a + b — c 
2t + to 
b + c — a 
2t + to 


= 0, 


= 0. 


(4.18) 


The situation is similar to the KdV/cKdV case. In the absence of time dependent terms (formally 
to —t oo), Eq. (14.181) reduces to a diagonal system. 

We take initial values of a, b , c to be 


a(v, 0) = 


l ?>o! b ^ =1 r c ^ = °- 


(4.19) 


Corresponding to the initial condition (14.191) (see also Fig. [TJ) , the Whitham system (14.18[) for 
the classical BO equation (formally to —> oo) admits an exact rarefaction wave solution in terms of 
the self-similar variable £ = r)/t. Itis6=l/2,c = 0 and from 


(2a-Oa'(O = 0 


(4.20) 
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Figure 7: Initial values (14.1911 for Riemann variables a, b and c. 


we find a = a(£) = £/2 (cf. 0); we also note that in j7(|, stationary solutions of Whitham type 
systems for the BO-Burgers equation were investigated. Like the cKdV equation there is a self 
similar system of equations which can represent the solution of the cBO equation for a portion of 
the domain; this system is given by 

2 (a — ftj + a + b — c = 0, 

2b l { b ~^j + a + b ~ c = °> (4.21) 

2Cj + b + c — a = 0, 

where £ = r]/ (2 t + to). 

As in the cKdV case, finding an analytical solution of the general Whitham system (|4.18l) is more 
difficult because of its non-diagonal property. Therefore we follow the numerical approach developed 
for the cKdV case in order to understand the structure of DSWs in the cBO equation. 


4-3. Comparison between numerical simulations of Whitham modulation equations and direct 
numerical simulations of BO/cBO 

In order to numerically solve the Whitham modulation equations in terms of Riemann type 
variables (|4.18|l we need to first obtain the boundary conditions. For the BO equation the boundary 
conditions are fixed in time and can be read from Eq. (14.191) or Fig. [7] In the case of the cBO 
equation (as opposed to the cKdV equation) we can find the boundary conditions for the Riemann 
variables analytically. Indeed, neglecting the derivatives with respect to the spatial variable ry in the 
Whitham system (14.181) . we obtain a linear system whose solutions can be easily found. The exact 
solution for the boundary conditions on the left side with the initial conditions (14.191) are given by 


a_ 

c_ 


R 2 {t) R{t ) 

2 2 

R 2 (t ) R{t) 

2 2 

\ \m - 1 ] 


[!-*(*)]- 

[1 - R(t )}, 


1 

2 ’ 


(4.22) 


where R{t) is given by (12.201) . Similarly the boundary conditions on the right side corresponding to 
the initial conditions (14.191) are found to be 


a + =b + = ^, c+ = 0. (4.23) 

As in the cKdV case, all Riemann variables for the cBO equation at the boundaries decay in time; 
see Fig. |8] for plots of these variables for t < 30. 
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(a) (b) 

Figure 8: Evolution of Riemann variables for cBO case which are (a) given by (14.2211 at the left boundary, 
(b) given by (14.231) at the right boundary. 


In order to obtain the numerical solutions of the Whitham system (14.181) . as with the KdV/cKdV 
case we again use Shampine’s hyperbolic PDE solver with the same version of the Lax-Wendroff 
method. In the numerical computations, we use N = 2 14 points for the spatial domain [—30,30] 
with the time step being 0.9 times the spatial step. The results computed with the initial condition 
(14.191) and the boundary conditions given above are shown in Fig. [9] at time t = 7.5. We note that 
there appear to be derivative discontinuities at the leading and trailing edges of the DSW. These 
may be smoothed by keeping higher order terms, but doing so is outside the scope of this paper. 




rt r/ 

(a) (b) 

Figure 9: Riemann variables at t=7.5 which are found by numerical solutions of reduced Whitham system 
for BO eq. and exact Whitham system (13.241) for cBO eqs. (a) for BO eq., (b) for cBO eq. Here, we take 
to = 10. 


We can reconstruct asymptotic solutions of the DSWs for both BO and cBO at any time t 
from the numerical solutions of the Riemann variables using Eqs. (14.161) and (14.171) with 9 given by 
Eq. (13.271) . These reconstructed solutions are compared with direct numerical simulations in Fig.HOl 
Next we solve the BO and cBO equations (i.e. Eq. (12.191) with c = 0 and c / 0 respectively) 
numerically with the initial condition (12.61) regularized as in Eq. (13.351) . The numerical method used 
is completely analogous to KdV/cKdV, except that e 2 d vvv is replaced by e'Hd riri . In this computation, 
we use N = 2 14 spatial Fourier modes with the domain size L = 30, and choose the time step to be 
10” 4 . The regularization parameter in the initial condition (13.351) is chosen to be K = 10, and the 
parameters in the cBO equation (12.191) are taken to be to = 10 and e = 10 -3 / 2 . Direct numerical 
simulations of BO/cBO and solutions of the associated Whitham equations are compared in Fig. flOl 
at t = 7.5. In the reconstructions from the Whitham equations, we fix the arbitrary constant phase 
9q in Eq. (14.171) by adjusting the maxima of the Whitham reconstruction to agree with those from 
direct numerical simulations. We provide another perspective by including space-time plots of the 
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(a) (b) 


Figure 10: Numerical and asymptotic solutions of BO and cBO eqs. at t=7.5 with the initial data (IQIl . 
(a)for BO eq., (b) for c.BO eq. Here, we take to = 10 and e = 10~ 3 / 2 . 


direct numerical solutions of BO and cBO eqs.- see FigfTTl 

From Fig. |T0]it is clear that the results of direct numerical simulations and those of the Whitham 
equations are overall in very good agreement. While the humps in the leading edge of the DSW are 
captured nearly perfectly, we note that at the trailing edge there are discrepancies. This is due to 
the fact that as time increases the DSW humps move to the right and spread apart (unlike KdV). 
This phenomena is also observed in Fig[lT]for both BO and cBO eqs. For large time the trailing edge 
of the DSW becomes small and we expect higher order terms in the Whitham modulation equations 
will need to play a role. As with determining the arbitrary phase 6q, the higher order analysis is 
outside the scope of this paper. 
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Figure 11: Space-time plot of the direct numerical solutions between t = 0 and t = 20 (a)for BO eq., (b) for 
c.BO eq. Here, we take to = 10 and e = 10 -3 ^ 2 . 


For the BO equation, we observe that the amplitudes of the leading hump increase slowly in 
time; they are calculated to be 3.431, 3.586 and 3.648 for t = 10, 20 and 30, respectively. This 
suggests that as t —> oo the amplitude may asymptote to 4. The average speed of the leading hump 
of the DSW is approximately V avg = 0.831 at t = 7.5. This speed is close to the phase speed of 
the algebraic solitary wave solution of the BO equation [43| with an amplitude 4V avg = 3.34; this 
solitary wave is approximately represented by 




1 + 


4V 

V(rj-Vt) 


(4.24) 
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For the cBO equation, the average speed of the leading hump of the DSW at t = 7.5 is 
approximately V avg = 0.438. This speed is considerably smaller than in the BO case because 
the amplitude of the leading edge decreases in time. The trailing edge of the DSW looks similar to 
that of the BO equation but its amplitude also decreases in time (see Fig. USD- These features are 
similar to those observed for DSWs in the cKdV equation. 


5. Comparison with direct numerical simulations of KP/2DBO 


In this section we will solve the KP and 2DBO equations numerically using a modified version 
of Trefethen’s code (Program 27 in |44|) and compare with direct numerical simulations of the 
corresponding (1+1) cylindrical equations. First Eq. (12.41) is written in 2D Fourier space as 


u t + Cu+ -k x u 2 = 0. 


(5.1) 


Here u and u 2 are Fourier transforms of u(x : y,t ) and u 2 (x,y,t), respectively and the form of the 
operator C for the KP and 2DBO equations are 


- iXk 2 „ , 

C = -^ - ie 2 kl 

k x + A 


and 


C = 


iXk 2 „ 

:-^ -*esgn (k x )k x . 

k x + A 


(5.2) 


(5.3) 


We note that a regularization parameter A is added to the denominator C to prevent the singularity 
in (ED near k x = 0. Then by using integrating factor e tc , Eq. ED is written in the following 
equivalent form 


( e ' f “), 




+£,,2 _ 


u 2 = 0 . 


(5.4) 


Finally, a fourth order Runge-Kutta method is used for time integration of Eq. (15.41) . For Fourier 
spectral methods, the initial condition must be periodic and smooth. However, the parabolic front 
initial condition (12.51) does not satisfy these conditions. Therefore we use the following regularized 
initial condition instead of Eq. (12.51) in the numerical computations 


ui{x, y, 0) = - 


p tanh ( x + 


P(y, o) 


— /ztanh + lo + 


P(y, o) 


exp — tyi‘ 


iP 


2 y P 

Ly 

(5-5) 

Here the 2D computational domain is [— L X ,L X \ x [— L y ,L y \. This initial condition decays for large 
x ; consequently we have effectively imposed periodic boundary conditions. We choose the location 
of the backward front to be far from the forward front and the curvature of the backward front 
P(y , 0) to be much smaller than the curvature of the forward front P(y , 0); this minimizes the effect 
of the backward front. The parameter K regularizes the fronts (forward and backward), while the 
parameters m and p smooth the initial condition in the y-direction (see Fig. 1121) . 

In all numerical simulations we use large domain sizes L x and L y , and an initially parabolic front 
P(y, 0) = cy 2 with c = 0.1. We choose the regularization parameters K , l Q , m and p such that the 
relevant dynamics are locally equivalent between the exact initial condition (12.51) and the numerical 
initial condition (15.51) . 

The regularization parameter A that prevents the singularity at k x = 0 in Eq. m is chosen 
to be the complex number A = iXo , where Ao = 2.2204 x 10” 16 is the smallest floating number in 
computations which MATLAB allows 24]. Meanwhile, in Fourier space, we have ui (0, k y ) ^ 0 for all 


k y € R for the numerical initial condition (15.51) (this is also true for the exact initial condition (12.51) 1. 
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Figure 12: Contour plot of the numerical initial condition (15.51) for y, — —1 ,K = 10, P(y,0) = 0.1y 2 ,lo = 
40, P(y, 0) = O.Olj/ 2 , m = 0.05, p = 3, L x = 60 and L y = 20. 

Even though the term e tC near k x = 0 is a very small number, this nevertheless has an effect on 
the zero background of the solution which is the natural background of the solution: the magnitude 
of the background changes with time t. Therefore, to compare with the numerical solutions of 
1 + 1 dimensional cylindrical equations, we readjust the numerical solutions of the 2 + 1 dimensional 
equations to an appropriate background when we compare solutions. 

For the KP equation, we choose the spatial resolution to be 2 14 x 2 10 , the domain size to be 
L x = 60 and L y = 20, and the time step to be 10~ 4 . We use the parameters I\ = 10, Iq = 40, m = 
0.05 ,p = 3,P(y, 0) = O.Oly 2 and e 2 = 10 -3 . The numerical solutions of the cKdV equation with 
the ETDRK4 method and the KP-II equation with Trefethen’s code at y = 0 and y = ±1.25 are 
compared in Fig. [13] for t = 7.5. To provide another perspective, we give a contour plot of the 
numerical solution of the KP-II equation -see FigHT] also provided is an animation of the DSW 
propagation in KP-II between t = 0 and t = 8 (see [451] 4. 




Figure 13: . Comparison of numerical solution of cKdV equation and numerical solution of KPII equation 
for t=7.5 (a) at y = 0, (b) at y = ±1.25. Here, we take to = 10 and e 2 = 10~ 3 . 

The solution of the corresponding 1 ± 1 dimensional cylindrical equation coincides with the 

solution of the 2±1 dimensional equation only at y = 0. For the comparison of results at cross 

sections different from y = 0, the solution of the 1±1 dimensional cylindrical equation must be 

shifted in the horizontal direction by a value which can be determined by the solution of the FS 

equation (12.171) for the desired y-cross section and time t. Specifically we recall y = x + P(y,t)/ 2, 
~ 2 

where P(y, t) = with A = 1, c = 0.1. Therefore, for the comparison of solutions of cKdV and 

KP-II equations, shifting is performed by x = y — 0.03125 for y = ±1.25 and t = 7.5. 

DSWs in KP-II equation are observed in both FiglTTland the animation given in [iH ]. 
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Figure 14: Contour plot of the numerical solution of the KP-II equation at t = 7.5 with computation 
parameters for the KP-II eq. given in the text. 


For the 2DBO equation, we choose the spatial resolution to be 2 14 x 2 10 , the domain size to 


We use the parameters K = 10, 


be L x = 50 and L y = 20, and the time step to be 10 -4 . 
lo = 30, m = 0.05, p = 3, c = 0.1 and e = 10 -3 / 2 . 

The numerical solution of the cBO equation with ETDRK4 method and the 2DBO equation 
with Trefethen’s code at y = 0 and y = ±1.25 are compared in Fig. [15] for t = 7.5. Similarly in 
KP-II, contour plot of the numerical solution of the 2DBO equation at t = 7.5 given in FielTOland 
propagation of DSWs in 2DBO between t = 0 and t = 8 obtained by the numerical solution is 
achieved (see [46| for the animation). 

Here shifting is also performed by x = rj — 0.03125 for y = ±1.25 and t = 7.5; this is because the 
evolution of the front shape as given by the FS equation (12.171) are identical between KP and 2DBO. 
DSWs in 2DBO equation are observed in both FigfTbland the animation given in |46|. The spreading 
behavior of DSW humps of BO type equations that we mentioned section 4.3 is also observed for 
the time evolution of 2DBO eq. in the animation pbi ]. 




x x 

(a) (b) 

Figure 15: . Comparison of numerical solution of cBO equation and numerical solution of 2DBO equation 
for t=7.5 (a) at y = 0, (b) at y = ±1.25. Here, we take to = 10 and e = 10 -3 ^ 2 . 

In both comparisons, there is excellent agreement between solutions of the 2±1 dimensional 
equations and the 1±1 dimensional cylindrical equations at the chosen y-cross sections and time. 
This supports Whitham’s asymptotic method and the 1±1 dimensional cylindrical equations as 
accurate reductions of the 2 ± 1 dimensional equations.Indeed in terms of computational time, 
the numerical solutions of 2±1 dimensional equations, 1±1 dimensional equations and Whitham 
systems are on the order of days, hours and minutes, respectively. Hence, in addition to enhanced 
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understanding of the underlying DSWs, the asymptotic method has enormous advantages in terms 
of computational time. 

We note that the level of agreement between 2+1 and 1+1 numerical solutions at subsequent 
times depends on the y-domain size L y in the 2 + 1 dimensional equations. If the y-boundaries are 
relatively far away from y = 0, then initially the DSW does not hit the y-boundaries. However 
after a finite amount of time, the DSW reaches the y-boundaries and gets reflected. To delay this 
reflection effect, L y and correspondingly the number of grid points in y need to increase. 



X 


Figure 16: Contour plot of the numerical solution of the 2DBO equation at t = 7.5 with computation 
parameters for the 2DBO eq. given in the text. 

Another aspect of the numerical solutions of 2+1 dimensional equations to note is the evolution 
in the rr-direction. Since we use the regularized initial condition (15.51) . as time evolves oscillations 
in x exist on the left tail of the solutions. After some time, these oscillations affect the localized 
DSW behavior of the 2+1 dimensional equation and thus agreement with the solution of the 1+1 
dimensional cylindrical equation. To avoid the effect of these oscillations, Iq, L x and correspondingly 
the number of grid points in x need to increase for better agreement. 

6. Conclusion 

In this paper we consider DSW behavior and exact reductions across a parabolic front of 2 + 1 
dimensional KP and 2DBO equations to the 1 + 1 dimensional cylindrical KdV (cKdV) and the 1 + 1 
dimensional cylindrical BO (cBO) equations. We derive their associated modulation equations and 
write these equations in suitable Riemann coordinates. We solve the resulting Whitham systems 
numerically and compare these results with direct numerical simulations of the 1 + 1 equations. 
Apart from an unimportant phase and a small discontinuity in front of the DSW, the results are in 
excellent agreement. The discontinuity can be accounted for by higher order terms; but this detail 
is outside the scope of the present paper. We also compare the Whitham theory of cKdV with KdV 
and cBO with BO; we find that while the amplitudes of the DSW structures of KdV and BO remain 
0(1), the DSWs of cKdV and cBO decay slowly in time. We compare the DSW behavior across 
the parabolic front of the 2 + 1 systems to their 1 + 1 counterparts; after accounting for an small 
mean term; excellent agreement is obtained. We conclude that the cKdV/cBO equations are able to 
accurately describe DSW behavior along a ‘flattening’ parabolic front to the KP/2DBO equations. 

7. Acknowledgements 

This research was partially supported by the U.S. Air Force Office of Scientific Research, under 
grant FA9550-12-1-0207 and by the NSF under grant DMS-1310200 and by the Scientific and 
Technological Research Council of Turkey (TUBITAK) under grant 1059B-19140044. 


22 















References 


[1] J. Lighthill, Waves in Fluids, Cambridge University Press , Cambridge, UK, 1978. 

[2] N.F. Smyth, P.E. Holloway, Hydraulic Jump and Undular Bore Formation on a Shelf Break, J. 
Phys. Oceanogr. 18 (1988) 947-962. 

[3] P.M. Bazin, La Propagation Des Ondes, Mem. Pres. Acad. Sci., Paris, 19, 495, (1850). 

[4] R.J. Taylor, D.R. Baker, H. Ikezi, Observation of Collisionless Electrostatic Shocks, Phys. Rev. 
Lett. 24(5) (1970) 206-209. 

[5] M.A. Hoefer, M.J. Ablowitz, I. Coddington, E.A. Cornell, P. Engels, V. Schweikhard, Dispersive 
and Classical Shock Waves in Bose-Einstein Condensates and Gas Dynamics, Phys. Rev. A. 74 
(2006) 023623. 

[6] M.A. Hoefer, M.J. Ablowitz and P. Engels, Piston dispersive shock wave problem, Phys. Rev. 
Lett., 100, (2008) 084504. 

[7] W. Wan, S. Jia, J.W. Fleischer, Dispersive Superfluid-like Shock Waves in Nonlinear Optics, 
Nat. Phys. 3 (2007) 46-51. 

[8] C. Conti, A. Fratalocchi, M. Peccianti, G. Ruocco, S. Trillo, Observation of a gradient catastrophe 
generating solitons, Phys. Rev. Lett. 102 (2009) 083902. 

[9] J. Fatome, C. Finot, G. Millot, A. Armaroli, S. Trillo, Observation of Optical Undular Bores in 
Multiple Four-Wave Mixing, Phys. Rev. X. 4 (2014) 021022. 

[10] G.B. Whitham, Non-linear Dispersive Waves, Proc. R. Soc. A. 283 (1965) 238-261. 

[11] A.V. Gurevich, L.P. Pitaevskii, Nonstationary Structure of a collisionless Shock Wave, Sov. 
Phys. JETP-USSR. 38(2) (1974) 291-297. 

[12] P. Lax, C. Levermore, The Small Dispersion Limit of The Korteweg-De Vries Equation 1, 
Commun. Pure Appl. Math. 36(3) (1983) 253-290. 

[13] A.V. Gurevich, L.P. Pitaevskii, Averaged Description of Waves in the Korteweg- de 
Vries-Burgers Equation, Zh. Eksp. Teor. Fiz. 93 (1987) 871-880. 

[14] G.A. EL, R.H. J. Grimshaw, A.M. Kamchatnov, Evolution of Solitary Waves and Undular Bores 
in Shallow-Water Flows Over a Gradual Slope With Bottom Friction, J. Fluid Mech. 585 (2007) 
213-244. 

[15] Y. Matsuno, Nonlinear Modulation of Periodic Waves in the Small Dispersion Limit of the 
Benjamin-Ono Equation, Phys. Rev. E. 58(6) (1998) 7934-7939. 

[16] Y. Matsuno, V.S. Shchesnovich, A.M. Kamchatnov, R.A. Kraenkel, Whitham Method For the 
Benjamin-Ono-Burgers Equation And Dispersive Shocks, Phys. Rev. E. 75 (2007) 016307. 

[17] B.B. Kadomtsev, V.I. Petviashvili, On the Stability of Solitary Waves in Weakly Dispersing 
Media, Sov. Phys. Dokl. 15 (1970) 539. 

[18] M.J. Ablowitz, H. Segur, Long Internal Waves in fluids of great depth, Stud. Appl. Math. 62 
(1980) 249-262. 

[19] M.J. Ablowitz, P.A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering, 
Cambridge University Press , Cambridge, UK, 1991. 


23 


[20] F. Calogero, A. Degasperis, Solution by the Spectral Transform Method of a Nonlinear 
Evolution Equation Including as a Special Case the Cylindrical KdV Equation, Lett. Nuovo 
Cim. 23 (1978) 150-154. 

[21] R.S. Johnson, On the Inverse Scattering Transform, the Cylindrical Kortewg-de Vries Equation 
and Similarity Solutions, Phys. Lett. 72A(3) (1979) 197-199. 

[22] F. Kako, N. Yajima, Interaction of Ion-Acoustic Solitons in Two-Dimensional Space, J. Phys. 
Soc. Japan. 49 (1980) 2063-71. 

[23] X.P. Wang, M.J. Ablowitz, H. Segur, Wave Collapse and Instability of Solitary Waves of a 
Generalized Kadomtsev-Petviashvili Equation, Physica D 1 78, (1994) 241-265. 

[24] C. Klein, C. Sparber, P. Markovich, Numerical Study of Oscillatory Regimes in the 
Kadomtsev-Petviashvili Equation, J. Nonlinear Sci. 17(5) (2007) 429-470. 

[25] M.J. Ablowitz, W.C. Curtis, On the Evolution of Perturbations to Solutions of the 
Kadomtsev-Petviashvili Equation Using the Benney-Luke Equation, J. Phys. A. 44 (2011) 
195202. 

[26] C.Y. Kao, Y. Kodama, Numerical Study of the KP Equation for Non-Periodic Waves, Math. 
Comput. Simul. 82(7) (2012) 1185-1218. 

[27] A.T. Cates, D.G. Crighton, Nonlinear Diffraction and Caustic Formation, Proc. R. Soc. Lond. 
A. 430 (1990), 69-88. 

[28] P.N. Sionoid, A.T. Cates, The Generalized Burgers and Zabolatskaya-Khokhlov Equations: 
Transformations, Exact Solutions and Qualitative Properties, Proc. R. Soc. Lond. A. 447 (1994) 
243-270. 

[29] E.A. Zabolotskaya, R.V. Khokhlov, Quasi-plane Waves in Nonlinear Acoustics of Confined 
Beams, Sov. Phys. Acoust. 15 (1969) 35. 

[30] C. Lin, E. Reissner, H.S. Tsien, On Two-Dimensional Non-steady Motion of a Slender Body in 
a Compressible Fluid, J. Math. Phys. 27 (1948) 220. 

[31] J.C. Luke, A Perturbation Method for Nonlinear Dispersive Wave Problems, Proc. R. Soc. A. 
292 (1966) 403-412. 

[32] M.J. Ablowitz, D.J. Benney, Evolution of Multi-Phase Modes For Nonlinear Dispersive Waves, 
Stud. Appl. Math. 49(3) (1970) 225. 

[33] M.J. Ablowitz, H. Segur, On the Evolution of Packets of Water Waves, J. Fluid Mech. 92 (1979) 
691-715. 

[34] E. Hopf, The Partial Differential Equation ut + uu x = p xx , Comm. Pure Appl. Math. 3 (3) 
(1950) 201-230. 

[35] P.F. Byrd, M.D. Friedman, Handbook of Elliptic Integrals and Physicists (2ed.), 
Springer-Verlag , Berlin, 1971. 

[361 M.J. Ablowitz, D.E. Baldwin, Dispersive Shock Wave Interactions and Asymptotics, Phy. Rev. 
E. 87 (2013) 022906. 

[371 L.F. Shampine, Solving Hyperbolic PDEs in MATLAB, Appl. Numer. Anal. Comput. Math. 
2(3) (2005) 346-358. 

[38] B. Engquist, P. Lotstedt, B. Sjogreen, Nonlinear Fillters For Efficient Shock Computation, 
Math. Comp. 52 (1989) 509-537. 


24 


[39] M.J. Ablowitz, D.E. Baldwin, M.A. Hoefer, Soliton Generation and Multiple Phases in 
Dispersive Shock and Rarefaction Wave Interaction, Phys. Rev. E. 80(1) (2009) 016603. 

[40] S.M. Cox, P.C. Matthews, Exponential Time Differencing For Stiff Systems, J. Comput. Phys. 
176 (2002) 430-455. 

[41] A.K. Kassam, L.N. Trefethen, Fourth-order time-stepping For Stiff PDEs, SIAM J. Sci. 
Comput. 26(4) (2005) 1214-1233. 

[42] T.B. Benjamin, Internal Waves of Permanent Form in Fluids of Great Depth, J. Fluid Mech. 
29(3) (1967) 559-592. 

[43] H. Ono, Algebraic Solitary Waves in Stratified Fluids, J. Phys. Soc. Jpn. 39 (1975) 1082. 

[44] L.N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, USA, 2000. 

[45] Dispersive Shock Wave Propagation in KP-II Equation Between t = 0 and t = 8, 
youtu.be/AE xAQHRSjvE . 

[46] Dispersive Shock Wave Propagation in 2DBO Equation Between t = 0 and t = 8, 
youtu.be/aXU NYKFlkeO, 


25 



